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Abstract 

For present study, setting Strouhal Number (St) as control parameter, numerical simulations for 
flow past oscillating NACA-0012 airfoil at 10^ Reynolds Numbers (Re) are performed. Temporal 
profiles of unsteady forces; lift and thrust, and their spectral analysis clearly indicate the solution 
to be a period-1 attractor for low Strouhal numbers. This study reveals that aerodynamic forces 
produced by plunging airfoil are independent of initial kinematic conditions of airfoil that proves the 
existence of limit cycle. Frequencies present in the oscillating lift force are composed of fundamental 
{fs), even and odd harmonics (3/*) at higher Strouhal numbers. Using numerical simulations, 
shedding frequencies (fs) were observed to be nearly equal to the excitation frequencies in all the 
cases. Unsteady lift force generated due to the plunging airfoil is modeled by modified van der Pol 
oscillator. Using method of multiple scales and spectral analysis of steady-state CFD solutions, 
frequencies and damping terms in the van der Pol oscillator model are estimated. We prove the 
applicability of this model to all planar motions of airfoil; plunging, pitching and flapping. An 
important aspect of currently-proposed model is capturing the time-averaged value of aerodynamic 
lift coefficient. 
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1 Introduction 


To propose efficient and better designs for small swimming and flying unmanned vehicles, un¬ 
derstanding of the unsteady mechanisms to generate lift and thrust forces at very low Reynolds 
numbers (Re) is of key importance. Being a highly nonlinear system, fluid flowing over these 
vehicles carries great complexities. Since Knoller [10] and Betz [4] revealed the mechanisms for 
production of thrust due to vortex pattern behind oscillating and pitching airfoils, respectively, 
this field has attracted the attention of many researchers around the globe. Recently, due to the 
advent of biomimicking flying (micro air vehicles) and swimming robots (underwater vehicles), 
many research efforts in this direction are being presented. Ho et al. [7], Triantafyllou et al. [27], 
Wang [30], Lehmann [11] and Shyy et al. [20] have provided detailed reviews regarding progress in 
the field of unsteady aerodynamics for flapping flight of insects, birds and robots. Study of mecha¬ 
nism for generation of unsteady forces is still to be understood completely due to a wide spectrum of 
parameters that are involved. Both the experimental and currently available numerical techniques 
require costly resources involving a large amount in terms of time and money. Considering this fact, 
researchers have also focussed towards development of the reduced order models. These models 
are based upon the reduced number of states for a dynamical system. These models can be built 
by either phenomenological or direct model-reduction approaches. For phenomenological modeling 
technique, the behavior of a parameter from the response of a physical system is modeled by a set 
of ordinary differential equations. Self-excited oscillators are great examples of this type covering 
the range from electric circuits to the aero-elastic phenomena. On the contrary, a direct approach 
requires computation of coherent structures in the flow to develop a reduced-order model [22]. 
Proper-orthogonal decomposition (POD) based reduced-order models are classical examples in this 
category. It forms the reduced-order basis functions from the flow snapshots capturing optimal 
energy of a dynamical system. The governing equations, such as the Navier-Stokes equations, are 
projected onto these reduced basis to form a reduced-order model thereby reducing the system from 
millions of degrees of freedom to order of ten. These models have been successfully employed for 
flow control, design, optimization, and uncertainty quantification [1]. In the phenomena-based ap¬ 
proach for the development of reduced order models, work by Skop et al. [23] is of vital importance. 
They introduced the Skop-Griffin parameter in order to propose a modified van der Pol oscillator 
model to predict the lift of a cylinder coupled with a linear model for mathematical representation 
of structure’s motion. This model is a center of focus in the study of vortex-induced vibrations. In 
Ref. [24], they extended their model from rigid to elastic cylinders. Identifying the cubic nonlinear¬ 
ity in lift force and quadratic relation between lift and drag forces for a circular cylinder, Nayfeh et 
al. [16,18] used method of multiple scales [14,15,17] and proposed the first-order accurate reduced 
order self-excited oscillator models for steady and transient parts of the aerodynamic forces. Follow¬ 
ing the higher order spectral analysis technique [5,6] for the identification of nonlinear parameters, 
Qin [19] developed forced van der Pol oscillator models for the rotational, inline and transversal 
oscillations of circular cylinders considering primary resonance with soft and hard excitation cases 
at Re = 10^. Janajreh et al. [8] extended this study for rotational oscillations of circular cylinders 
at Re = 10^. Marzouk et al. [12] presented a second order accurate model for the steady-state lift 
and drag forces for stationary circular cylinders at different Re values. Keeping eccentricity as the 
control parameter and using method of harmonic balance, Akhtar et al. [2] performed the numerical 
simulations for flow over the elliptical structures and approximated the transient and steady state 
behavior of the lift and drag forces using a combined van der Pol-Duffing oscillator model. These 
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afore-mentioned studies help enhance our understanding regarding the nonlinear mechanism for 
production of the unsteady aerodynamic forces by the bluff bodies. To extend this technique for 
the aerodynamic bodies, Ellenrieder [28] adopted the methodology proposed by Shop et al. [25] 
to present following oscillator model for time-dependent lift produced by the airfoils undergoing 
forced plunging motion; 

Cl - cosG{Cl^^ - AQ^)Cl + oos^Cl = uj^Fh ( 1 ) 

where ojg is the fundamental vortex shedding frequency, Cl^ is the maximum amplitude of the 
time-dependent lift force coefficient, h is an instantaneous position of airfoil while plunging and, G 
and F are the constants determined from the empirical formulation proposed by Shop et al. [25]. 
In the current study, we first show the existence of limit-cycle behavior in the response of an airfoil 
performing forced oscillatory motion. Here, the response of this nonlinear system is determined 
in terms of lift and thrust forces. We consider different sets of initial kinematic conditions of the 
oscillating airfoil. Identifying the nonlinearities in the response, we present a modified forced van der 
Pol oscillator model for lift-force coefficient. The analytical solution of this nonlinear mathematical 
model is derived using method of multiple scales; a powerful perturbation technique. The strength 
of this technique to solve the perturbation problems lies in the fact that both the small and large 
values of system’s states can be handled by the involvement of slow and fast time-scales. This 
reduced-order model not only captures the temporal details of the the lift force but also it predicts 
its spectral composition accurately. The purpose of this study is to provide a computational tool 
for estimating the response of this physical system quickly. To determine the parameters of the 
presented reduced-order model, we employ numerical data obtained from CFD (computational 
fluid dynamics) simulations using ANSYS Fluent [3]. We perform these CFD simulations for a 
range of St values at Re = 10^. These simulations are thoroughly validated and used here as 
the first step towards development of the reduced-order models. The manuscript is organized as 
follows. Section 2 provides necessary details on our CFD simulation methodology for flow past the 
oscillating airfoils along with the validation studies. Using different sets of the kinematic initial 
conditions and states of this nonlinear system in section 3, we show the limit-cycle existence in 
the response of the oscillation airfoils. Mathematical formulation and derivation of a reduced- 
order model based on a modified forced van der Pol oscillator equation is presented in section 4. 
Applicability of this model for lift produced by plunging, pitching and flapping airfoils is shown 
in section 5. This model is capable of accurately predicting the nonlinear behavior of unsteady 
Cl oi oscillating airfoils. To model its time-averaged value Cl-, we introduce a different type of 
quadratic nonlinearity to the self-excited van der Pol oscillator model and describe its details in 
section 6. Section 7 explains the overall behavior of the linear and nonlinear damping terms with 
respect to increasing control parameter; the Strouhal number (defined in section 2). We also show 
the results of predictive settings for the present model to compute unsteady Cl and its spectra for 
the intermediate values of Strouhal number where the model parameters were calculated through 
interpolation of the data-set from CFD solutions. In sections, we conclude and summarize our 
present work. 

2 Numerical Methodology 

For the present study, we simulate flow past an oscillating NACA-0012 airfoils by solving two- 
dimensional incompressible Navier-Stokes equations using ANSYS-Fluent; a finite-volume based 
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commercial software. The Navier-Stokes equations in their integral form can be written as; 

/^ pcpv.dA - T^Vcp.dA- I S^dV = 0 

V V 

where p is the density of fluid, v represents the velocity vector, A shows the surface-area vector, 
r<^ is the diffusion term, Vcf) denotes the gradient term and shows the source term. Since the 
present case employs moving mesh technique [3], the source term is zero. The temporal term is 
approximated using the first-order implicit method. Second-order upwind scheme is employed for 
numerically approximating the convective term. Assuming incompressible flow presently, velocity 
and pressure terms are coupled through the pressure implicit with splitting of operator (PISO) 
algorithm. To avoid effect of disturbances on boundaries, radius of the 0-type domain is kept 
at 25c, where c is the chord-length of airfoil, as shown in Fig. 1. Flow domain is meshed using 
unstructured triangular cells. We use high grid resolution near the airfoil surface to resolve the 
boundary layer and to capture the wake characteristics downstream of the airfoil. Airfoil surface 
is resolved using 400 nodes. Dynamic meshing techniques; spring analogy and remeshing, are 
employed in the vicinity of moving airfoil for grid transition that allows the adjustment of the 
grid in accordance with the instantaneous position of airfoil. Numerical solutions of flow-helds are 



Figure 1: Schematic of Geometry and Fluid-Domain 

highly dependent on the suitability and accuracy of boundary conditions. In Fluent, motion of an 
object may be defined by a user-defined function ( UDF) which is a computer code written in the 
C-language environment coupled with the Fluent-Macros. Forced motion of an airfoil can be of 
three types; plunging, pitching and flapping (combination of plunging and pitching). The schematic 
of these motion are shown in Fig. 2. Pithing and plunging motion can be modeled as; 

Plunging : h(t) = ho cos(27r/f -|- cph) 

Pitching : a{t) = Uo sin(27r/f) 

where ao is the maximum pitching amplitude and (ph denotes the phase-angle. Dirichlet conditions 
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Figure 2: Airfoil Kinematics (Solid Lines: Upstroke, Dotted Lines; Downstroke) (For clarity, only 
mean and the amplitude positions are shown here) 

are employed on the inlet boundary and the pressure outlet condition is used for the outflow 
boundary. At this boundary, static pressure is specified. For incompressible flows, pressure on the 
boundary is determined by taking the average of specified values on the cell faces and computed 
values of the static pressure on the corresponding cell-centers. All other flow variables are computed 
by extrapolation of the computed values from the inner domain. The Reynolds number, defined 
as Uooc/h', is calculated using appropriate values of the kinematic viscosity u, keeping all other 
parameters equal to unity. Strouhal number (St) is defined as St = 2fho/Uoo where is the 
maximum amplitude, / is the excitation frequency for oscillating airfoil in Hertz, and Uoo is the 
free-stream fluid velocity. It is considered as the primary governing parameter for investigating 
unsteady behavior of oscillating airfoils, ho depicts the wake-width behind a plunging airfoil. For 
pitching motion, total vertical excursion traversed by the trailing-edge of airfoil approximates the 
wake-width. St is varied by changing plunging and pitching amplitudes while incoming reference 
flow velocity and oscillation frequency are kept fixed. In the present study, numerical simulations 
are initialized by the uniform inlet velocity boundary conditions. To perform these simulations, we 
employ 1.18440 x 10® number of cells in the whole domain and 2000 time-steps per oscillation cycle 
of the airfoil. Details for the grid-convergence, time-step refinement, and validation studies are 
available in Ref. [9]. Lift and thrust forces are computed by integrating pressure and shear stresses 
over the surface of airfoil. Their corresponding coefficients. Cl and Ct respectively are computed 
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as; 




L 


nPUoo C 


( 4 ) 


Ct = y- 


^pUoo C 

Using time-period of an oscillation cycle (r = 1//), corresponding time-averaged coefficients are 
calculated through following relation; 


1 

C= C(t)dt 

TJt 


( 5 ) 


3 Existence of Limit-Cycle 

Limit cycles can be used to model various nonlinear oscillatory systems from the real world. 
Limit-cycle exists if energy input to the dynamical system is balanced by an energy dissipation 
mechanism and the system’s response limits itself onto either a periodic, a quasi-periodic, or a 
chaotic attractor [13]. Hence this phenomenon represents a closed trajectory in the phase space. 
To characterize the behavior of a nonlinear system as a limit cycle, it needs to be tested for 
various initial conditions. Based upon a reduced-order model for Cl, von Ellenrieder [28] and von 
Ellenrieder et al. [29] presented the aerodynamic response of plunging airfoil as a limit cycle. To 
obtain numerical solution for governing ordinary differential equation, all the initial conditions were 
described in the form of and its time-derivatives. Use of the actual initial kinematic conditions 
seems more appropriate to investigate the limit-cycle behavior of any dynamical system. For current 
research-work, we consider a broader spectrum of the initial conditions through inclusion of initial 
position where the airfoil starts its motion from. There can be infinite starting positions between 
positive and negative amplitudes of an oscillating airfoil. In this study, initial conditions refer to the 
initial kinematic states of plunging motion for the numerical simulations. These include following 
scenarios; 

1. Does plunging airfoil undergo upstroke or downstroke first? 

2. Which position does it start its oscillation from? 

We choose three positions to start the airfoil’s oscillatory motion. These are the top-most, mean 
(mid) and the bottom-most positions. Airfoil can perform upstroke initially starting from the 
bottom-most or mean position and it may undergo downstroke starting from mean or the top¬ 
most position, (fill = 0° gives initial downstroke starting from the maximum positive amplitude of 
oscillating airfoil. Using these conditions of initial stroke and position, four cases are required to 
be studied for one set of values of oscillation frequency and amplitude. To analyze this phenomena 
here, we consider four initial conditions described earlier in terms of the starting position and 
stroke of a plunging airfoil. We compute the response of this dynamical system in terms of the 
aerodynamic lift and thrust coefficients denoted as Cl and Ct, respectively. Using the temporal 
histories and phase maps of these dynamical states, we observe similarity in their nonlinear character 
for different initial conditions. In Fig. 3, Cl and Ct for four initial conditions of plunging airfoil 
are presented. We observe the same phase difference in their temporal histories as given in the 
forced motion of plunging airfoil. Another notable feature is the same amplitude of for all the 
initial conditions. Phase of is exactly equal to that of heave motion. But for Ct, phase depends 


6 







Figure 3: Comparison of Cl and Ct for different initial positions and strokes of plunging airfoil for 
St = 0.10 in (a) and (c); and for St = 0.30 in (b) and (d) 


upon magnitude of 4> only and not its sign (where we can say 270° = —90°). Here, we support 
the existence of limit-cycle using aerodynamics forces and their time-derivatives as the states of 
flow past a plunging airfoil through phase maps shown in Fig. 4 for St = 0.10 and 0.30. Periodic 
steady-state solutions for the aerodynamic force coefficients are presented for this purpose. Like 
bluff-body aerodynamics [2,12,18], the frequency of unsteady Ct is twice of that for the unsteady 
Cl- All the initial conditions lead us to the same period-n attractor thus proving independence of 
this system from initial kinematic states of oscillating airfoil for the range of parameters considered 
here. 

4 Reduced Order Modeling 

Young et al. [31] explained the existence of vortex-lock in phenomena for plunging airfoil for a 
range of flow and kinematic parameters. In such cases, natural vortex-shedding frequency comes 
out to be equal to forcing frequency of oscillating airfoil. Presence of higher-order harmonics in 
spectra of unsteady aerodynamic force coefficients for ffow over oscillating airfoils exhibits this 
phenomena as nonlinear. This is a well-known fact in vibrations of bluff-bodies like cylinders and 
cables as well [21,26]. In current study, spectra of the unsteady lift force are analyzed to identify 
the prominent frequency components and their corresponding harmonics present in the signal. For 
given range of St, (from 0.05 to 0.5), fundamental vortex shedding fg frequency appears to be equal 
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Figure 4: Phase maps using Cl and Ct their time-derivatives for different initial positions and 
strokes of plunging airfoil at St = 0.10 in (a) and (c); and for St = 0.30 in (b) and (d) 


to the forced plunging frequency (n/27r), equal to 0.5Hz of the airfoil. It indicates that these all 
are the cases related to the primary resonance. First even and odd harmonics show prominent 
peaks. As we increase St, even harmonic starts getting larger amplitude. Figure 5 shows three 
shedding cycles represented by unsteady Cl and its relevant amplitude-spectra. The amplitudes 
of energetic fundamental, first even and odd harmonics are also indicated being employed here for 
model-reduction. An asymmetric forced van der Pol oscillator model is used here to analytically 
represent the lift force produced by plunging NACA-0012 airfoil. Due to the presence of even 
harmonics in C^-spectra, we introduce a quadratic nonlinear term containing a multiple of ClCl- 
Choice of this term to model the quadratic nonlinearity is made due to a phase-difference of nearly 
7 r /2 (or its odd integral multiples) between the fundamental and first even harmonics. Assuming 
this system as weekly damped and having a soft excitation; 

Cl + sCl = ~ olClCl ~ iCl^Cl + F’oCos(Dt -|- P)] (6) 

Us is the vortex-shedding frequency, e is a book-keeping parameter while /r, a and 7 represent the 
linear, quadratic and cubic damping coefficients, respectively. All these parameters are positive 
real numbers. Although no vortex-shedding from NACA-0012 is observed at zero angle-of-attack 
at this Re, we consider this case related to primary resonance (D ~ Us) in a mathematical sense; 


Q. = u + ea 


( 7 ) 













Figure 5: (a) Unsteady Lift and (b) Amplitude-Spectra for St = 0.30 

Using Eq. 7 and Fo as the amplitude, excitation may be expressed as; 

E{t) = Focos[{i0s + ea)t+ T] (8) 

Assuming time scales are To = t, Ti = et and T 2 = e^t, we have; 

T'(t) = To cos[(a;sTo-I-crTi + r] (9) 

Using Chain rule for differentiation, we can write; 

— = Do + eTi -|- D 2 + ... (10) 

at 

(f 

— = D\ + 2eDoDi + e2(T\ + 2DoDi) + ... (11) 

Let second-order approximate solution for Eq. 6 be; 

C'L(t) = CLo(To,Ti,T2) + eC'Li(To,Ti,T2)+e2CL2(To,Ti,T2) (12) 

We expand all the terms in Eq. 6 as follows; 

Cl = (T^o + 2eDoDi e^D^i 2 ^DoD2)(Clo + ^Cli + ^^Cl2) 

= D\Clo + eD\CLi + e^D'^oCL2 + 2 €DoD,Clo + 2e^DoDiCLi 

-|- D'^iClo + 2e^D 0 D 2 CL 0 (13) 

oP'sCl = 0}^ s{Clo + (-Cli + ^^Cl2) 

= Oj'^sClo + €<^‘^sCli + sCl2 (14) 

^IJ‘Cl = Gfi{Do + tDi + D 2 + ■■■){Clo + ^Cli + ^Cl2) 

= enDoCio + f^DiCii + e^/^TiCio (15) 
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(16) 


—eaCiCL — ~^oi{Clo + fC*Li + c^Cl2)(-^o + ^Di + ^ Di){Clo + fCii + c^Cl2) 
= —eaCioiDoCLo) — ^^aCioiDoCLi) 

— e^aCLoiDiCio) — e^aC'Li(-DoC'Lo) 


-eiCLCt^ = -e^[CLo + eCLi + e^CL2f{Do + eDi + ^Di){Clo + ^Cli + 

= -eiCL\{D.CLo) - ehCL\{D.CLo) 

- e^CL^iDiCLo) - 2e^-fCLoCLi{DoCLo) (17) 

eFoCos{uJsTo + (tTi +r) remains in its present form. Equating like powers of e on both sides of 
Eq. 6, we get; 

0(e°): D^oClo+u:sClo = ^ (18) 


0{€^) : D'^oCli + ujsCli = —2DoDiClo + I^DoClo — oiClo{DoClo) — iCl^o{DoClo) 

+ Fo cosiuigTo + crTi + E) (19) 


0{^) : D\Cl2+^sCl2 = -2DoDiCli - D\Clo - 2D,D2Clo + - aCUD.CLi) 

— aCLoiDiCio) — oiCLiiDoCLo) 

- iCL^iDoCLi) - 7 Cl\{DiClo) - 2^CloCli{DoClo) ( 20 ) 


Solution of Eq. 18 is; 

Cloit) = Ao{Ti,T2)cOs[u}sTo + 9o{Ti,T2)] 

To solve Eq. 19, all the terms on right hand side need expansion. 

^ O'^Clo _ N , o, . ^ ^ , a ^ 

^ dT dT — Qj-^ sinyujgTo H“ ^o) “h qj-^ cosi^ujsTo ~\~ Oo) 


( 21 ) 


( 22 ) 


= -fiuJsAo sm{usTo + 6o) 


(23) 


= ^auJsAo'^ sm(2u}sTo + 20o) 


(24) 


-iCl. 


2dCLo _ 1,„ , /I 3 ■ ' ~ " ' 1 


dT ^ sin(a;5To + 9o) + ^'yujgAo^ sm{3uJsTo + 39o) 


(25) 


Fo cos{uJsTo + crTi + E) = Fo cos(a;sro + 9o + uTi + r - 9o) (26) 

= Fo cos{uJsTo + 9o) cos{aTi + T - 9o) - Fo sin(cjsro + 9o) sin(crri + E - 6*o) (27) 
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Modulation equation may be obtained by combining multiples of sin^LOgTo + 9o 


-Fo sin(cjrir0o) = 0 

ul 1 4 


( 28 ) 


dAo nAo 'jAo^ Fo 

- 7 ^ 7 ^ =-smlcrTi + F + Po) 

dTi 2 8 2a;^ ^ ^ 

By combining multiples of cos(a;sro + 60 ), second modulation equation may be constructed. 

do 

2oJsAo-^ = Fo cos(itTi + r - 0 o) 
oil 

To make the system autonomous, we assume; 

r] = crTi + r - 6*0 
9o = crTi + T - rj 

dOo dr] 


So, Eq. 30 may be written as; 


dr] , ^° ] I 


Thus amplitude and phase are governed by; 


■ ]iAo 'jAo^ Fo . ^ ^ ^ ^ ^ 

^° = -2 -+ ' 


(29) 


(30) 


(31) 


(32) 


?7 = cr + 


-cos{r]) 


2 iijgAo 

To get the steady-state motion, time-derivative of both amplitude and phase needs to be zero. 

]iAo ^Ao^ F 


(33) 


= -—sin((TTi -I- r -I- 0o 


8 2 ,oj 


(34) 


(^Ao = 7r^cos(r/) 
ZiOs 


Squaring and adding both Eq. 34 and 35, we get; 


Fo , ]rAo aAo ,2 , /i 2 2 
= (-^ +Ao a 


^ 2 

Identifying the modulating terms, we get following governing expression for CLi{t); 

d^Cr^ A A^ 

—^ -h uIs'^Cli = auJs-^sm{2uJsTo + 29o) + 'yuJs-^sm(3u}sTo + 39o) 

dTo 2 4 


(35) 


(36) 


(37) 
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To get particular solution, we can write; 


Clip = K sin(2a;sT’o + 2do) + M sin(3a;sro + Wo) 


dC 


= 2ujsK cos(2uJsTo + 29o) + WgM cosISw^Tq + 39o 

(jTq 


(38) 

(39) 


d‘^C 


= —Aujs'^K sin(2cjsTo + 29o) — 9 l0s‘^M sin(3wsTo + 39o) 

dTo 


(40) 


putting these values of Clip, and ^ into Eq. 37 and equating like terms, we get K = 

and M = —Hence, Ciiit) comes out to be; 


cxA ^ 'yA ^ 

CLipit) = - „ ° sin( 2 a; 5 t + 29o) - sin(3a;5t + 39o) 


(41) 


Thus complete solution for Ciit) is; 

aA. 


TT 27^0 




Ciit) = Ao cos{u)st + 9o) - e— -cos( 2 a; 5 t + 29o + w) - t— — cos{3i0st + 39o + -) + ••• (42) 


Gcu 


32cUo 


For autonomous dynamical system in complex notation, we have; 

CL{t) = + g-d^+r-*?)] 

_ T i(2nt-^2T—2p-^^) _j_ —^(2fit+2r—2?7 +^)i 

12 w^ ^ ^ 

_ 7^0 r i(3nf+3r-3r)+g^) , -j(3f2t+3r-3r)+^)i 

64a;, ^ ^ ^ 


(43) 


Now, we describe the procedure for identification of linear and nonlinear parameters. Performing 
the spectral analysis of unsteady lift profiles, values for fundamental vortex shedding frequency u, its 
corresponding amplitude Ao and amplitudes yl 2 and A^ of first even and odd harmonics, respectively 
can be measured. Comparing the assumed approximate solution of van der Pol oscillator model in 
Eq. 12 with that from multiple scales method in Eq. 42, nonlinear damping parameters; a and 7 , 
can be calculated as; 


a = 


6 a;, H 2 


(44) 


External detuning parameter is; 


7 = 


32a;s2l3 


(45) 


H — a; 


a = 


(46) 
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Setting phase F of the excitation, t] can be measured from the phase (p of component in 

Fourier Transform of lift signal; 


v = r-cp[CLm 


Forcing amplitude Fo is; 


liOs^oO' 

cos(r/) 


Linear damping fi can be calculated as; 


7 ^ 0 ^ FoSin(?7) 

4 AoUs 


(47) 


(48) 


(49) 


5 RESULTS & DISCUSSION 

When airfoil starts oscillation in a fluid, a reverse von Karman vortex street is observed in the 
wake which leads to generation of time-varying aerodynamic forces. At low St, these aerodynamic 
forces are periodic in nature which may be composed of several harmonics of fundamental vortex 
shedding frequency. We compare the solution of van der Pol oscillator model with those from CFD, 
both in temporal and spectral domains. The proposed model not only captures temporal prohles 
accurately but also strong harmonics, in their spectra. We present the validity of the proposed 
reduced-order model for three types of forced motions of airfoil; plunging, pitching and flapping. 


5.1 Plunging Airfoil 

Plunging airfoil performs oscillatory motion along vertical direction. It portraits a single-degree 
of freedom system. Figure 6a shows comparison of aerodynamic lift force from CFD simulations 
and the proposed model in time domain. Solution for the ordinary differential equation (ODE) of 
proposed reduced-order model is obtained by numerical integration using Runge-Kutta method of 
order-4. It is clear from these figures that current model fulhlls the requisites effectively. To capture 
the inherent nonlinear characteristics of lift force, the results of the proposed model should carry 
the same spectral components of the signal. It can be analyzed by comparing the spectra of signals 
from CFD and reduced-order model. Ci-spectra of CFD and van der Pol oscillator model are 
shown in Fig. 6b. It shows that the proposed model not only captures the fundamental frequency 
with accurate amplitude but also the other harmonic components are well predicted by this model. 
To develop a database for linear and nonlinear damping parameters, force amplitude, forced and 
vortex shedding frequencies. Table 1 presents sample dynamic parameters for a range of St. 

5.2 Pitching Airfoil 

Although pitching motion of airfoil is a different degree-of-freedom but its response in terms of 
aerodynamic forces resembles those of plunging airfoil. Table 2 presents sample values for spectral 
and ROM parameters for different values of St. In Fig. 7, we show comparison of CFD solution 
with that of van der Pol oscillator model. 
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Figure 6: Comparison of CFD solution and the proposed model for Plunging Airfoil at St = 0.45 
(a) Time Histories (Solid Lines: CFD and Circles: ROM) (b) Amplitnde-Spectrum (Solid Lines: 
CFD and Dashed Lines: ROM) 


Table 1: Spectral and ROM Parameters for Lift of Plunging Airfoil 


Parameters/St 

0.15 

0.30 

0.45 

fs 

0.4982 

0.5003 

0.4998 

Ao 

1.908 

6.051 

10.3 

A 2 

0.001039 

0.01653 

0.06951 

A^ 

0.08816 

0.3482 

0.2933 

V 

37.0033° 

48.2008° 

53.0429° 


1.1382 

1.4396 

0.7144 

a 

0.0054 

0.0085 

0.0123 

7 

1.2714 

0.1580 

0.0270 

Fo 

0.1760 

0.1513 

0.0870 


5.3 Flapping Airfoil 

As described earlier, flapping motion is a combination of plunging and pitching motions. Prac¬ 
tically, swimming and flying species or robots employ synchronized pitching and plunging while 
flapping their wings. Looking at the similar nature of response, we model lift force of flapping 
wings using same oscillator model. Table 3 shows sample values of spectral and ROM parameters 
for flapping airfoil while Fig. 8 shows comparison of CFD and ROM solutions. CFD simulations 
presented here were carried ont for 0.05 < /lo < 0.50, = 10°, / = 0.5Hz and Re = 10^. R also 

proves the suitability of this model for this complex phenomenon. 
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Table 2: Spectral and ROM Parameters for Lift of Pitching Airfoil 


Parameters/St 

0.10 

0.20 

0.30 

fs 

0.5002 

0.5001 

0.5003 

Ao 

0.8367 

2.049 

3.576 

A2 

0.0005242 

0.0008243 

0.002885 

A3 

0.0136 

0.1203 

0.3186 

V 

-0.4177° 

6.376° 

10.8025° 


0.4098 

1.4757 

2.2212 

a 

0.0141 

0.0037 

0.0043 

7 

2.3351 

1.4059 

0.7008 

Fo 

0.0072 

0.0081 

0.2209 




Figure 7: Comparison of CFD solution and the proposed model for Pitching Airfoil at St = 0.30 
(a) Time Histories (Solid Lines: CFD and Circles: ROM) (b) Amplitude-Spectrum (Solid Lines: 
CFD and Dashed Lines: ROM) 

6 Improved Model 

Like usual unsteady signals, we may decompose Cl into two components; time-averaged value 
Cl and fluctuating component Cl 


Cl = Cl + Cl (50) 

Although motion of airfoil during plunging, pitching and flapping is symmetrical about its mean 
position, yet there exists a non-zero time-averaged value of [32], Presence of non-zero time- 
averaged value of a signal along with an even harmonic shows quadratic nonlinearity in the system 
[17]. In section 4, we attempt to model the quadratic nonlinearity using a term ClCl- It is chosen 
due to a phase-difference of tt/2 (or its integral multiple) between fundamental and first even 
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Table 3: Spectral and ROM Parameters for Lift of Flapping Airfoil 


Parameters/St 

0.05 

0.20 

0.35 

fs 

0.5017 

0.5003 

0.4999 

Ao 

0.6721 

2.093 

5.24 

A 2 

0.02461 

0.0288 

0.04316 

A 3 

0.03737 

0.07439 

0.1265 

ri 

-31.6684° 

60.7861° 

60.1451° 


1.4077 

0.8865 

0.6364 

a 

1.0304 

0.1240 

0.0456 

7 

12.4164 

0.8161 

0.0924 

Fo 

0.0467 

0.0542 

0.0450 




Figure 8: Comparison of CFD solution and the proposed model for Pitching Airfoil at St = 0.20 
(a) Time Histories (Solid Lines: CFD and Circles: ROM) (b) Amplitude-Spectrum (Solid Lines: 
CFD and Dashed Lines: ROM) 


harmonic in the unsteady Cl signal. This model can capture the unsteady details effectively but 
it does not predict Cl that is an important feature of aerodynamics in case of oscillating airfoils. 
Considering lesser magnitudes of forcing function amplitude Fq, a possible singularity sX r] = tt!2 
and to take care of Cl, we propose another version of the modified van der Pol oscillator’s ordinary 
differential equation to model the lift of oscillating airfoil. 

Cl + oj\Cl = eiiiCL - aCl - ^Cl^Cl] (51) 
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Solving this system using the method of multiple scales, following second-order solution is obtained; 


Ciit) = Ao cos{u}t) + e[ 


aAo"^ aAo^ 


2uP‘ 6a;2 


cos( 2 ti;t -|- 2/3o 


aAo" 
~ 32a; 

Modulation equations in this case are; 


cos(3a;t -|- 3/3o)] 


(52) 



io = l^Ao 

0 

(53) 


■ _ 117^^0"^ 5a^Ao^ 

~ 8 a; 8 a; 256a; 12a;3 

(54) 

Linear and nonlinear damping parameters are identified as; 



8 a; A 3 

A. 

(55) 


3 cu ^ A 2 

(56) 


4 /i 

V 

(57) 

From Eq. 52, Cl is; 

^ 2 a/r 

Cl = 9 

70;^ 

(58) 


To prove the accuracy of this model, we compare its results with those from CFD in Fig 9. CFD 
results for St=0.20 and 0.40 were used here to extract the ROM parameters. The strength of this 
model lies in its capability to model time-averaged value of the lift force. This model, too, behaves 
well in time as well as spectral domains. We also show a database for ROM parameters for different 
St values in Table 4. Comparing damping parameters with those in Table 1 for lift of plunging 
airfoil, we observe almost equal magnitudes. It is a manifestation for capability of this improved 
reduced-order model to predict actual unsteady characteristics in lift force signal. 


7 Predictive Settings 

To summarize the results of the proposed model and its suitability for responses of various 
kinematics of airfoil, Fig. 10 shows the variations of damping coefficients as functions of St. Positive 
values of ;U, a and 7 indicate the presence of limit cycle. Despite a difference in level of their 
magnitudes in case of three different kinematics of airfoil, similar trend for variation can be observed 
for quadratic (a) and cubic ( 7 ) damping-coefficients. This similar nature of nonlinear dynamic 
parameters proves the effectiveness of the present reduced-order model. Magnitudes of nonlinear 
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(b) St=0.40 

Figure 9: Comparison of solutions from CFD and the improved ROM for Plunging Airfoil (a) Time 
Histories (Solid Lines: CFD and Circles: ROM) (b) Amplitude-Spectrum (Solid Lines: CFD and 
Dashed Lines: ROM) 
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Table 4: Parameters for Lift of Plunging Airfoil in Case of Improved ROM 


Parameters/St 

0.15 

0.30 

0.45 


1.1613 

1.4462 

0.7157 

a 

0.0085 

0.0134 

0.0194 

7 

1.2760 

0.1580 

0.0270 

Cl (ROM) 

0.0046 

0.0316 

0.1083 

Cl (CFD) 

0.0063 

0.0352 

0.1190 





□ Heaving Airfoil 
V Pitching Airfoil 
O Flapping Airfoil 


Figure 10: Plots of Dynamic Parameters for Varying St 


damping parameters decrease due to increasing amplitude Ao for fundamental harmonic as can be 
seen in Eq. 44 and 45. Linear damping /r for of pitching airfoil comes out to be a linear function 
of St while it shows quite similar nonlinear behavior for plunging and flapping motions. Now, we 
examine the performance of our reduced-order models in predictive settings for St = 0.18. For this 
purpose, the selection of this St value is justified from Fig. 10 where we observe higher gradients in 
the trends of both linear and nonlinear damping parameters. We find out the linear, quadratic and 
cubic damping values using data points presented in Fig. 10 through cubic interpolation scheme. 
We compute the solution of the proposed van der Pol oscillator model by numerical integration 
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and compare its performance with the CFD results. Figure 11 shows that the reduced-order model 
well satisfies the requirement not only in time-domain but also in the spectral domain even at the 
point around highest variation in the available data. To quantitatively test the model, we calculate 
the percentage errors in the amplitudes of unsteady Cl from CFD and predicted values of reduced 
order model, and the fundamental frequency in the spectra of both. These values come out to be 
7.91% and zero, respectively. While comparing both the spectra, we see small deviation in the 
levels of higher harmonics. Cl comes out to be 0.009 from CFD and the numerical solution of the 
model in predictive settings gives 0.006 which are quite close to each other. Errors and deviations 
from actual CFD results may be removed by using more data points as samples in the interpolation 
scheme. 




Figure 11; Comparison of CFD solution and the proposed model for Pitching Airfoil at St = 0.20 
(a) Time Histories (Triangles: CFD and Solid Line: ROM) (b) Amplitude-Spectrum (Dashed Line: 
CFD and Solid Lines: ROM) 


8 Conclusions 

In the present paper, we identify the existence of a limit-cycle behavior through aerodynamic 
forces produced by oscillating airfoils. It is done through various initial kinematic states that has 
greater physical significance. Noting the presence of even and odd harmonics in C^-spectra for 
plunging, pitching and flapping airfoils, we propose a phenomenological reduced-order model by 
employing an asymmetric forced van der Pol oscillator model. Parameters for this models are 
calculated through results of CFD simulations. By numerically integrating the governing ROM 
ODE, we present comparison of its solutions with those of CED. Results are quite promising for 
both temporal and spectral domains. Real strength of this proposed model lies in its applicability 
to aerodynamic lift forces of plunging, pitching and flapping airfoils that are altogether different 
degrees-of-freedom of the relevant structure. Similarity of nonlinear behavior in all of these phe¬ 
nomenon are also proven by plotting nonlinear damping coefficients versus Strouhal number. This 
model helps measure the aerodynamic forces of oscillating streamlined body quickly and accurately 
without involvement of costly experimental equipment or time-consuming complex CED simula- 
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tions. This model is limited in a sense that it cannot predict Cl. To overcome this deficiency and, 
considering the lower amplitudes of forcing functions and a probable singularity condition in forced 
van der Pol oscillator model, we propose another model that carries Cl^ term to model quadratic 
nonlinearity with no forcing function. This model captures Cl quite effectively. Magnitudes of 
linear and nonlinear dampings in Table 3 and 4 are quite close to each other. 
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